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1.  INTRODUCTION 


For  the  first  quarterly  report^  we  showed,  in  admittedly  brief  terms, 
how  the  code  INTERN  is  constructed,  and  the  physical  basis  for  its  content. 

In  order  to  demonstrate  the  viability  of  the  code  some  very  preliminary  cal¬ 
culations  for  both  a  towing  tank  and  the  real  oceanic  environment  were  performed. 

During  the  quarter  being  reported  upon  here  INTERN  has  undergone  sub¬ 
stantial  evolution.  Also,  attempts  at  validation  of  the  code  have  improved 
our  understanding  of  the  physical  and  numerical  basis  of  the  code  a  great 
deal.  This  validation  process  has  included  comparison  with  experiment,  which 
is  the  ultimate  test  of  any  theoretical  analysis.  So  far,  this  comparison 
against  experiment  has  been  very  successful  as  the  contents  of  this  report 
will  reveal. 

In  an  effort  to  generalize  the  code,  provision  has  been  made  for  the 
inclusion  of  a  sail  for  a  submerged  body.  This  work  is  also  described  below. 

In  previous  work  it  has  been  possible  to  assume  left-right  symmetry,  but  if 
one  wishes  to  include  shear  currents  or  swirl,  such  a  symmetry  condition  is 
not  applicable.  Initially  INTERN  was  designed  under  the  assumption  of  such 
symmetry.  This  is  no  longer  necessary,  and  INTERN  is  now  able  to  treat  problems 
involving  shear  and  swirl. 

Also  implemented  in  the  code  is  the  capability  to  make  plots  of  the 
crosstrack  surface  current.  Such  plots  are  useful  to  obtain  a  pictorial  idea 
of  the  surface  pattern  that  is  generated  by  the  passage  of  a  body  below  the 
surface. 

Below  we  present  details  on  the  above  subjects  as  well  as  a  discussion 
of  the  sample  calculations.  A  report  from  our  subcontractor  Flow  Research 
on  their  far  field  code  is  also  included. 
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Two  problems  in  particular  have  us  concerned  about  the  proposed  far  field 

calculations:  one  problem  is  that  of  storage,  which  may  likely  be  overcome;  the 

other  problem  has  to  do  with  the  modal  approach  itself.  The  modal  approach, 

it  seems  to  us,  has  virtue  only  if  the  number  of  modes  that  contribute  is  not 

large,  TJiat  is,  convergence  is  obt  tined  with  only  a  relatively  small  number  of  modes 

ci  atributing.  Present  codes  generally  include  twenty  modes.  Calculations  by 

(2) 

Selwyn  indicate  that  a  significant  fraction  of  the  energy  can  reside  in  the 
higher  modes,  so  that  it  is  difficult  to  have  confidence  that  twenty  modes  (or 
even  substantially  more)  result  in  convergence.  Included  in  this  report  are 
two  possible  suggestions  for  avoiding  the  difficulties  with  the  modal  approach 
presently  used  for  the  "analytical"  calculations. 
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2.  CODE  GENERALIZATIONS 


In  this  section  we  describe  some  of  the  changes  and  alterations  which 
have  been  made  to  INTERN  during  this  quarter.  Naturally  only  the  major  changes 
are  discussed  here,  since  at  this  stage  the  code  is  still  undergoing  rapid 
development. 


2.1  INCLUSION  OF  SAIL 

Bodies  which  are  designed  for  subsurface  travel  usually  have  a  sail  or 
conning  tower  placed  near  the  top  forward  position  of  the  body.  While  not  making 
up  a  large  fraction  of  the  body  volume,  it  is  possible  that  such  an  object  could 
contribute  to  or  change  the  internal  wave  pattern  generated  by  the  body. 

The  body  has  so  far  been  represented  in  INTERN  by  a  simple  source-sink 
combination  (Rankine  Ovoid).  At  a  later  stage  more  complicated  distributions 
will  be  utilized  to  more  precisely  represent  realistic  body  shapes.  The  sail^ 
can  be  modelled  in  potential  flow  by  means  of  a  combination  line  source  and 
line  sink  on  the  y  =  0  plane  (see  the  first  quarterly  report  for  the  definition 
of  our  coordinate  system)  parallel  to  the  z-axis. 

Let  m  be  the  source  (sink)  strength  per  unit  length  in  the  z-direction; 

then 


u(x)  ■  <x_x  ) 

-+ 

where  x'  is  the  position  vector  to  the  line  element  dz  .  In  our  case  x 
so  that 
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These  numbers  can  be  checked  by  locating  the  stagnation  point  for  the  cylindrical 
streamline,  ip  =  0  ,  around  the  line  source  and  sink,  and  are  assumed  to  apply 
for  a  body  1006  cm  at  maximum  diameter,  and  9915  cm  in  length. 

By  performing  the  indicated  integrals,  one  finds  for  the  velocity  field. 
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While  the  above  formulae  appear  to  be  a  trifle  messy,  they  are  really  very  simple 
to  include  in  the  subroutine  for  the  potential  flow.  No  parameter  studies  have 
been  made  to  determine  just  how  much  influence  the  sail  does  have,  but  a  sample 
calculation  has  been  made  which  indicates  that  the  sail  does  not  contribute 
significantly  to  the  internal  wave  generation.  When  time  permits  a  careful  com¬ 
parison  calculation  can  be  made.  Presently  the  sail  is  either  included  or 
excluded  in  the  calculation  by  means  of  a  control  number  ISAIL.  Naturally,  when 
the  distances  from  the  sail  are  such  that  the  sail  potential  flow  velocities  are 
very  small,  this  part  of  the  overall  potential  flow  is  bypassed  in  the  interest 
of  computational  speed. 
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The  result  that  the  sail  does  not  contribute  greatly  to  the  internal 
wave  generation  can  be  understood  physically  by  the  fact  that  the  fluid  tends 
to  go  around...the  nail  ra.th.er  than  be.  lifted,  up  over  it.  There  is  some  lifting 
near  the  top  of  the  sail,  but  only  for  a  short  axial  distance.  Also  perhaps 
one  ought  to  use  the  theory  developed  for  joined  bodies,  since  it  may  not  be 
possible  to  make  a  clear  separation  of  the  potential  flow  as  we  have  done  so 
far.  This  whole  problem  of  a  precise  potential  flow  will  be  discussed  more 
fully  at  a  later  date.  In  any  case  a  careful  parameter  study  needs  to  be  done 
before  a  firm  conclusion  about  the  influence  of  a  sail  is  drawn. 

2.2  ASYMMETRICAL  PROBLEMS 

The  original  design  for  INTERN  included  the  assumption  of  left-right 
symmetry  for  the  problems  of  interest.  This  assumption  is  unduly  restrictive 
and  a  considerable  effort  was  undertaken  to  eliminate  this  particular  restriction. 
A  thorough  debugging  was  essential  to  make  certain  that  the  code  still  gave 
correct  results.  It  is  now  a  relatively  simple  matter  to  run  problems  with  or 
without  left-right  symmetry.  A  control  number  LSY  /  0  is  used  to  designate 
asymmetrical  problems. 

WAKE  SWIRL 

Generally*  a  propeller  imports  a  certain  amount  of  angular  momentum  into 
the  wake  of  a  body.  An  efficient  propeller  will  minimize  the  amount  of  energy 
that  goes  into  angular  motion.  This  angular  or  rotational  motion  in  a  wake 
is  usually  referred  to  as  swirl.  It  would  be  convenient  to  have  a  model  for 
swirl  which  could  be  readily  used  in  a  research  tool  such  as  INTERN.  If  one 
neglects  the  axial  variations  of  swirl,  so  that  the  motion  is  strictly  two  di¬ 
mensional,  then  the  model  presented  in  Ref.  (3)  for  a  simple  type  of  eddy  is 
probably  an  adequate  one  for  our  purposes. 

Let  us  assume  that  the  swirl  then  depends  at  any  downstream  time  on  only 
the  radial  distance  from  the  track.  If  uq  is  the  orbital  velocity,  then  the  swirl 
motion  is  described  by 

Ar  -rz/4at 
u  =  - t  e 
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where  O  will  be  taken  as  a  turbulent  viscosity,  A  is  a  constant,  and  t  the  downtrack 
time  from  the  propeller.  Both  the  constant  A  and  the  initial  distribution  can  be 
obtained  from  the  assumption  for  the  initial  energy  associated  with  the  swirl 
motion,  which  should  be  something  less  than  10%  of  the  initial  wake  energy  (potential 
energy) . 

No  calculations  have  yet  been  made  using  the  above  model  for  the  swirling 
motion  behind  a  propeller  with  INTERN.  Such  calculations  would  be  very  easy  to 
perform  if  ARPA  so  desires.  It  is  felt  that  it  would  be  best  to  first  validate 
the  code  without  too  many  effects  included  so  that  any  problems  encountered  can 
be  identified.  When  confidence  is  assured  that  the  code  gives  correct  answers, 
then  we  can  readily  proceed  to  include  these  extra  effects.  For  example,  we 
should  be  satisfied  that  the  turbulence  model  is  satisfactory  before  straining 
its  credibility  in  the  presenci  of  current  shear  and  swirl.  Also,  it  is  important 
that  the  far  field  be  computed  satisfactorily  before  we  attempt  to  refine  the 
physical  content  of  the  code.  Another  point  is  that  the  presence  of  current 
shear  requires  the  nonlinear  terms  or  special  treatment  this  will  effect  the 
far  field  calculations. 

2.3  RESTART  CAPABILITY 

It  is  highly  desirable  that.it  be  possible  .to  restart  INTERN  in  case 
one  wishes,  on  any  particular  calculation, . to  compute  further  downstream.  This 
bookkeeping  task  has  been  accomplished,  and  it  is  possible  now  to  perform  such 
a  restart.  Many  features  of  this  nature  must  be  included  to  make  INTERN  a 

-I 

convenient  user's  code. 
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3.  COMPARISON  WITH  EXPERIMENT 


As  part  of  the  validation  process  for  INTERN  it  is  important  to  compare 
the  code  results  not  only  with  other  types  of  calculations,  but  with  experiment 
where  possible.  The  turbulence  model  had  been  compared  with  Naudascher's  data 
in  the  previous  quarterly,  but  during  this  last  quarter,  an  effort  was  made  to 
check  INTERN  against  experimental  results  taken  at  APL. 

Basically,  the  APL  experiment  consists  of  a  tank  containing  stratified 
water  through  which  a  self-propelled  12:1  slender  body  is  run.  The  body,  except 
for  the  absence  of  a  sail  and  trim  planes,  is  shaped  realistically.  In  our  cal¬ 
culation  we  utilized  a  Rankine  Ovoid  with  the  same  maximum  diameter  and  the  same 
length  as  the  APL  model.  The  specifications  for  the  comparison  are  shown  in 
Fig.  (1)  along  with  a  picture  of  the  computational  mesh.  The  stratification 
was  assumed  linear  in  the  calculation,  whereas  the  experiment  showed  some 
deviation,  particularly  very  near  the  top  and  the  bottom. 

The  calculation  was  run  twice,  once  without  the  wake,  i.e.,  utilizing 
body  generated  waves  only,  and  again  allowing  for  both  body  and  wake  generated 
i  waves.  As  in  our  comparison  with  Naudascher,  the  wake  is  initiated  four  diameters 
behind  the  body. 

Figures  (2)  and  (3)  show  computer  plots  of  the  density  variation  for  the 
case  of  body  waves  only.  Due  to  the  variable  time-step  in  the  wake  calculation 
we  do  not  have  comparable  plots  when  the  wake  is  present.  This  plot  could  easily 
be  added  for  the  wake  also.  For  the  actual  comparison  with  experiment,  the 
numbers  were  printed  out  and  hand  plotted  on  the  same  graph  paper  used  for  the 
experiment . 

Figure  (4)  demonstrates  the  comparison  that  was  obtained  between  theory 
and  experiment  which  appears  to  be  quite  good.  Although  we  understand  there  is 
more  complete  data  being  taken,  the  experimental  data  we  have  plotted  becomes 
much  less  reliable  beyond  60  diameters  downstream.  The  calculation  was  first 
run  using  data  assuming  the  same  turbulence  profile  that  was  used  to  compare 
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with  !•  ludascher '  s  data.  This  profile  made  the  calculation  become  out  of  phase 
with  ( xperiment  beyond  40  diameters  downstream.  By  using  a  value  qQ  based  on 
lu'l  -  .087  u  in  close  agreement  with  the  FRI  and  VPI  data  for  momentumless 

ii  o’ 

propeller  wakes  (Naudascher  used  a  jet),  the  phase  of  the  calculation  agrees 
very  well  with  the  APL  experiment. 

In  the  calculation  shown  in  Fig.  (4)  a  completely  mixed  wake  was  assumed, 
although  various  runs  using  widely  different  values  for  the  mixing  had  only  very 
small  effect.  Presumably  a  greater  effect  would  be  seen  downstream.  A  plot  of 
the  crosstrack  surface  current  downstream  of  the  body  is  shown  in  Fig.  (5).  This 
plot  is  made  by  connecting  values  of  the  surface  velocity  v^  in  the  y  direction 
at  many  axial  positions.  The  extreme  smoothness  of  the  plot  shows  that  the 
numerical  procedures  in  the  code  work  very  well  indeed.  Figures  (6)  and  (7)  show  some 
velocity  vector  plots  of  the  flow  field  in  the  tank.  We  observe  that  beyond 
100  diameters  downstream  there  are  severe  top  and  bottom  reflections  in  the  tank. 

There  are  several  points  that  should  be  made  about  this  comparison 
between  theory  and  experiment.  First  there  is  the  question  of  the  difference 
in  body  shapes  between  that  used  in  the  experiment  and  of  the  Rankine  Ovoid 
assumed  in  the  calculation.  This  could  account  for  the  difference  observed 
in  the  region  near  the  body.  Additionally,  while  the  experimental  angle-of- 
attack  is  limited  to  .004  radians  it  is  still  possible  that  such  an  angle-of- 
attack  could  help  explain  any  variance  between  theory  and  experiment.  Both 
of  these  effects  can  be  taken  account  of  in  INTERN.  Particularly,  the  angle- 
of-attack  question  should  be  addressed  for  very  practical  reasons.  We  have 
given  some  thought  to  this  question  and  are  confident  that  this  can  be  included 
in  the  primary  potential  flow.  Another  interesting  question  is  the  point  at 
which  the  wake  should  be  initiated.  If  this  point  is  chosen  poorly  then 
the  choice  will  have  an  adverse  effect  upon  the  computed  phases  of  the  internal 
wave  field.  Presently  INTERN  initializes  the  wake  at  four  diameters  -  which 
might  be  too  far.  This  point  should  be  investigated  more  thoroughly.  Another 
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point  in  regard  to  our  turbulence  model  is  the  fact  that  we  likely  must  allow 
for  the  suppression  of  the  vertical  turbulent  velocity  fluctuations  w' ,  i.e., 
our  model  for  turbulence  should  be  anisotropic.  Such  a  model  is  included  in 
WAKEVM,  and  can  be  incorporated  in  INTERN.  There  is  evidence  that  within  a 
Vaisala  period  the  turbulence  is  strongly  anisotropic. 

The  calculation  for  the  APL  tank  has  been  carried  farther  downstream  - 
to  210  diameters  (4.5  Vaisala  periods),  but  as  yet  we  have  no  data  to  compare 
it  to.  This  run  was  made  utilizing  realistic  density  profiles.  Hopefully 
we  will  be  able  to  make  comparisons  with  data  in  the  next  quarterly  report. 
Such  a  comparison  will  also  show  the  adequacy  of  our  present  turbulence  model 
farther  downstream. 

At  the  present  time  we  are  making  calculations  on  the  density  profiles 
sent  to  us  by  Dr.  Phil  Selwyn.  So  far  we  have  run  three  of  the  cases  and  have 
encountered  absolutely  no  difficulty  in  making  these  runs.  The  difficulty 
for  INTERN  (which  we  trust  is  very  temporary)  is  of  course  the  far  downstream 
or  far  field  calculations  which  we  discuss  in  more  detail  in  the  next  section. 
Since  INTERN  is  not  susceptible  to  the  same  kinds  of  problems  that  these  test 
problems  were  intended  to  show  we  do  not  plan,  in  the  interest  of  economy, 
to  run  all  the  test  problems. 
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Figure  2.  Computer  plot  of  density  variance  at  a  point  13.97  cm 
directly  above  body  centerline  (body  only) . 


PROBE  1 


Figure  4.  Comparison  of  theory  and  experiment  for  the  displacement  of  a  fluid  particle. 


X/D  -  4.85482E+01 


X/D  =  8.79183E+01 
Figure  6. 

Velocity  plots  for  APL  experiment 


X/D  -  1.01304E+01 


X/D  «=  1.53666E+02 


X/D  =  2.06028E+02 
Figure  7. 

Velocity  plots  for  APL  experiment 
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4.  FAR  FIELD  CALCULATIONS 


In  this  section  we  discuss  the  situation  regarding  the  far  field  cal¬ 
culations.  Below  in  Section  4.1  is  a  report  from  our  subcontractor  Flow  Research. 

As  mentioned  in  the  introduction,  the  far  field  approach  taken  by 
Flow  Research,  and  used  by  others  for  these  kinds  of  calculations,  namely  the 
modal  approach,  has  us  sufficiently  concerned  that  we  are  looking  at  other 
methods  for  far  field  calculations.  In  particular,  the  questions  of  completeness, 
convergence  and  computer  storage  are  the  main  items  of  interest.  In  Section  4.2 
we  present  an  alternative  analytical  approach  which  is  being  considered.  This 
approach  was  suggested  by  P.  Roache  at  SAI,  and  is  a  variation  upon  the  FRI  method. 

Another  method  being  considered  is  analogous  to  the  Region  I  solution 
in  INTERN.  This  method  has  the  advantage  of  being  simple,  and  of  being  accurate. 
Basically,  the  method  consists  of  turning  off  the  turbulence  and  the  potential 
flow  solutions  of  Region  II,  allowing  more  storage  capacity  and  a  great  increase 
in  speed.  Another  method  is  also  being  studied,  but  it  is  too  preliminary  to 
report  at  this  time. 

4.1  FRI  FAR  FIELD  CODE 

4.1.1  Summary  of  Progress 

During  this  report  period,  a  link  with  the  LBL  computing  facility  has 
been  established,  the  basic  elements  of  the  FRI  internal  wave  code  have 
been  developed  or  acquired  and  successfully  tested,  and  the  structure  for 
the  remaining  development  has  been  defined.  The  current  version  of  the  code 
calculates  the  internal  waves  generated  in  an  ocean  or  a  fluid  with  side- 
walls  by  a  moving  source-sink  disturbance.  It  has  been  partially  validated 
by  comparisons  with  computat-ions  based  on  the  work  of  Smithmeyer  and  Pao^ 
in  which  a  constant-N  fluid  in  assumed.  Validation  for  more  general  stratifi¬ 
cations  is  in  progress. 


18 


The  source-sink  calculation  is  based  on  the  following  expression: 


2  V  “  Sln(oT“) 

■’y,z)  tt  U  J  ^  2  *n  (_h,B  )(VZ,B  ) 


cospydB  (1) 


where  ij;  is  fluid  displacement  and  definitions  of  the  other  symbols  can  be 
found  in  Ref.  (5).  For  Eq.  (1)  and  in  general,  our  wave  code  utilizes 
Milder' s^  subroutines  for  computing  the  eigenmodes  of  the  equation 


subject  to  the  boundary  conditions 


4»(0)  =  <p  (-H)  =  0  . 


,  (7) 


Selwyn  s  improvements  to  Milder 's  routine,  ZMODES,  have  been  included. 

In  our  use  of  the  discrete  Fourier  transform  to  evaluate  inverse  integral 

Fourier  transforms  (see  Eq.  (1)  and  Ref.  (8)),  we  use  procedures  for  cal- 
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culating  real  sine  and  cosine  transforms  described  by  Cooley  1  to  take 
advantage  of  the  usual  circumstance  that  the  integrand  is  a  real,  odd  or 
even  function,  such  as  in  Eq.  (1).  The  speed  of  our  FFT  algorithm  is  thus 
increased  by  a  factor  of  4. 

4.1.2  Structure  for  Remaining  Development 

Early  in  the  development  of  our  code,  we  decided  to  lead  up  to  the  full 
complexity  of  obtaining  the  flow  downstream  of  a  plane  of  initial  data  by 
solving  a  sequence  of  simpler  problems.  This  approach  facilitates  the  job 
of  validation,  since  each  block  of  code,  starting  with  the  basic  eigenmode 
and  Fourier  transform  routines,  can  be  independently  tested. 

Having  solved  the  source-sink  problem,  the  next  step  in  our  code  develop¬ 
ment  plan  is  to  solve  the  initial  plane  problem  given  an  initial  distribution 
of  fluid  displacement  (the  wake  collapse  problem).  We  will  first  choose  a 
displacement  distribution  which  can  be  analytical] y  Fourier  transformed,  such  as 
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Hy,?.) 


where  L  and  1,  are  length  scales.  The  resulting  expression  for  displacement 

2.  y 

downstream  of  the  initial  plane  is 
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dz 


where  ¥  is  the  transform  of  the  initial  data.  This  version  of  the  code  will 
be  validated  by  comparing  its  results  to  those  of  Dugan^  ^  and  Young^  \  We 
then  will  allow  for  the  specification  of  an  initial  grid  of  displacement  data 
requiring  numerical  transformation.  Finally,  we  will  solve  the  full  initial 
plane  problem  by  permitting  the  specification  of  an  initial  grid  of  the  per¬ 
turbations  w,  v,  p,  and  p.  The  first  of  the  above  steps  will  require  a 
z-integration  procedure,  the  second  step,  in  addition,  will  require  a  procedure 
to  numerically  Fourier  transform  the  initial  data,  and  the  last  step  will 
require  a  generalization  of  these  procedures. 

An  important  option  to  be  added  to  the  code  as  time  permits  is  the  use 

of  a  stationary  phase  approximation  when  the.  distance  downstream  becomes  so 

large  that  the  discrete  transform  becomes  inefficient  because  of  wake 

spreading  and  increasingly  high  oscillation  frequencies  in  the  lateral 

spectrum.  In  fact,  the  stationary  phase  approximation  becomes  more  accurate 

as  the  discrete  transform  becomes  more  inefficient.  For  the  source-sink 

til 

problem  the  stationary  phase  result  for  the  integral  of  the  n  mode  is 
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=  0  . 


where  B  is  a  solution  to 
ns 


_d_ 

dB 


p / co s  9 

VYn 


+  sin 


For  additional  details  on  the  current  status  and  further  development  of  the 
code,  see  Section  4,1.3  of  this  report, 

A  final  note  is  that  many  of  the  details  of  the  interfacing  with  the 
near  wake  calculation  of  SAI  have  been  worked  out.  Most  importantly,  it 
has  been  decided  to  keep  the  SAI  and  FRI  codes  essentially  separate.  The 
FRI  code  will  obtain  the  initial  data  that  it  requires  from  a  tape  generated 
by  the  SaI  code.  This  modular  approach  avoids  the  complication  of  using 
overlays,  thereby  increasing  the  flexibility  and  transferability  of  the 
combined  SAI-FRI  programs. 

4.1.3  Details  of  Code  Structure 

Figure  8  is  a  block  diagram  of  the  FRI  internal  wave  code.  Existing 
modules  are  denoted  by  solid  lines,  whereas  dotted  lines  are  used  to  indicate 
future  development.  Each  arrow  corresponds  to  the  calling  of  a  subroutine. 
The  modules  are  described  below: 

EXISTING  ROUTINES 

-  controls  I/O  and  program  flow 

-  secs  up  the  array  defining  N(z) 

-  computes  cubic  spline  coefficients 

-  calculates  interpolated  values  using  coefficients 
computed  by  SPLINES 

-  produces  printer  plots 

-  controls  calculation  of  dispersion  relations  and 
eigenfunctions 

obtains  convergence  of  eigenvalue  -  eigenfunction 
iteration  procedure 


FRIWAC 

TCLINE 

SPLINES 

SPLINE 

GRAPH 

ZMODES 

EIGENVL 
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EIGENFN 

FINT 

CMODE 

SSINV 

SSNTRP 

SIX 

FIVE 

EXSET 

FFFT2 

FUTURE  ROUTINES 
IPLANE 
ZING 

IPINV 

SSSPV 

IPSPV 

CSPV 


-  uses  finite  differences  to  compute  trial  eigenfunctions 

-  does  third  order  interpolations  in  conjunction  with 
GRAPH 

-  controls  wave  field  calculation 

-  computes  contribution  to  wavefield  from  source-sink 
disturbance  using  discrete  Fourier  transforms 

-  interpolates  from  eigenmode  wave  number  grid  to  wave 
number  grid  required  for  discrete  Fourier  transforms 

-  takes  into  account  that  the  integrand  (function  of 
wave  number)  to  be  transformed  is  even 

-  takes  into  account  that  the  integrand  is  real 

-  sets  up  exponentials  for  use  of  FFT  algorithm 

-  uses  FFT  algorithm  to  compute  discrete  Fourier  transform 


-  will  Fourier  transform  a  grid  of  data  given  at  initial  plane 

-  will  compute  z-integrals  from  transformed  initial  plane 
data 

-  will  compute  contribution  to  wave  field  from  initial  plane 
data  using  discrete  Fourier  transforms 

-  will  compute  contributions  to  wave  field  from  source-sink 
disturbance  using  method  of  stationary  phase 

-  will  compute  contributions  to  wave  field  from  initial 
plane  data  using  method  of  stationary  phase 

-  will  calculate  stationary  phase  points 

-  will  provide  plots  of  the  wave  field 


PLTFLD 


IPLANE 


TCLINE 


ZMODES 


CMODE 


PLTFLD 
- T 


Figure  8.  Block  diagram  of  Flow  Research  internal  wave  code. 

Solid  lines  indicate  current  status,  dotted  lines 
denote  future  development. 
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FAR  FIELD  CALCULATIONS  OF  SURFACE  CURRENT 


In  this  section  we  describe  an  alternative  approach  to  the  far  field 
analytical  calculations  being  undertaken  by  FRI.  A  method  is  described  for  cal¬ 
culating  the  far  field  surface  current  produced  by  a  submarine  and  its  wake.  Initial 
data  are  provided  at  a  vertical  plane  in  the  wake  from  either  finite-difference 
calculations  including  turbulence  modeling,  or  from  an  analytical  potential-flow 
representation  of  the  submarine.  The  method  is  based  on  numerical  Fourier 
transforms  in  the  cross-track  direction,  numerical  Laplace  transforms  in  the 
track  direction,  and  a  finite-difference  solution  in  the  vertical  direction. 
Contrasted  to  a  previously  suggested  approach,  the  present  method  involves  no 
eigenfunction  problem,  and  requires  computer  storage  of  only  two-dimensional 
rather  than  three-dimensional  arrays. 

4.2.1  Introduction 

Reference  (1)  is  a  progress  report  on  the  problem  of  calculating  surface 
waves  from  a  submarine  or  other  moving  submerged  body,  including  turbulence 
in  the  near  wake.  The  flow  around  the  body  and  the  near  wake  are  to  be 
obtained  by  finite-difference  calculations,  as  described  in  Ref.  (1),  or  may 
alternately  be  described  with  less  accuracy  by  potential  flow  methods,  as 
in  Ref.  (12).  (However,  the  present  method  would  not  be  restricted  to  a  simple 
point  source  as  in  Ref.  (12),  but  could  use  arbitrary  distributions  of  sources, 
sinks,  and  circulation.) 

The  present  report  is  concerned  only  with  a  method  of  calculating 
the  far  field,  in  which  the  flow  is  assumed  inviscid.  The  theoretical  develop¬ 
ment  and  notation  is  that  of  Ref.  (1),  but  an  entirely  different  approach  is 
used  for  the  solution  of  the  linear  problem  of  the  far  wake.  The  body  moves 
along  the  positive  x-axis  (track  direction),  the  y-axis  is  in  the  cross-track 
direction,  and  the  z-axis  is  vertical,  with  z  =  0  at  the  surface  and  z  =  -H  at 
the  ocean  bottom.  The  far  field  (Region  III)  begins  at  a  distance  x^  behind 
the  body.  The  velocity  components  of  the  fluid  in  (x,y,z)  are  (u,v,w),  respectively. 
The  most  difficult  part  of  the  solution  is  to  obtain  9w/3z  at  the  surface. 
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The  point  of  departure  of  the  present  method  from  that  of  Ref.  (1) 
is  Eq.  (8)  of  Ref.  (1),  repeated  here. 


dz 


+  B 
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2n2 

a  U 


-  1 


W=  4  W*(B,z)  +  |  (B,z)  -  f- R-({p->'  E  F(a, B,Z)  (2) 

Mo  a*' 


This  is  a  complex  ordinary  differential  equation  in  z  for  w(a,B>z)  where  w 
is  the  doubly-transformed  w-velocity  component;  w(x,y,z)  has  been  Fourier- 
transformed  in  the  y-direction,  with  the  y-coordinate  transforming  to  the  real 
frequency  domain  3,  and  Laplace-transformed  in  the  x-direction,  with  the 
x-coordinate  transforming  to  the  complex  frequency  domain  a.  U  is  the  constant 
velocity  of  the  body,  and  N=N(z)  is  the  Brunt-Vaisala  frequency  given  by 


dp 

ZSl _ 2. 

p  dz 
o 


(3) 


where  g  is  the  acceleration  of  gravity  and  pQ(z)  is  the  density  of  the  fluid. 

No  restrictions  on  p  (z)  or  N(z)  are  required.  On  the  right-hand  side,  the 
*  *  *° 

terms  W  ,  V  ,  and  R  are  the  Fourier  transforms  in  y  of  the  initial  values  of 
(v.  v,  and  p  at  the  plane  x=x2>  i.e.  w(x2>y,z)  is  Fourier  transformed  to  W  (B,z), 
v(x2,y,z)  to  V  (B,z),  and  p(x2,y,z)  to  R  (B,z). 

In  Ref.  (1)  (pp.  24-25),  it  is  proposed  to  solve  this  problem  by 
expanding  w  in  eigenfunctions  which  must  be  solved  numerically.  The  question 
of  the  completeness  of  the  set  of  eigenfunctions  has  been  answered  only 
tentatively  (Ref.  (1),  pg.  26),  The  final  solution  still  involves  an  integration 
in  z  and  an  inverse  transform  in  B.  The  eigenfunction  solution  is  anticipated 
to  be  the  most  time  consuming  part  of  the  program,  and  will  require  500,000 
words  of  storage.  (The  CDC  7600  has  500,000  words  of  slow  core  storage  and 
70,000  words  of  fast  core  storage.) 

4.2.2  Suggested  Method 

Here,  we  propose  to  avoid  the  eigenfunction  problem  entirely.  For 
moderate-sized  programs,  we  could  do  all  calculations  in  fast  core,  with  some 
buffering  from  slow  core.  No  questions  arise  about  the  completeness  of  the 
eigenfunctions. 
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The  essence  of  the  suggested  method  is  to  solve  Eq.  (2)  numerically, 
and  to  collapse  a  storage  dimension  out  of  the  problem,  without  changing  the 
actual  dimensionality. 

Equation  (2)  is  an  o.d.e.  for  w(a,B,z).  The  right-hand  side  is 
three-dimensional,  F(a,3,z),  but  only  by  a  as  a  parameter.  The  complex 
quantities  W  ,  dV  /dz,  R  are  2-D  arrays  in  (B,z).  We  are  interested  only  in 
9w/9z  at  the  surface  value  of  z=0,  so  we  store  a  complex  2-D  array  9w/9z | z=Q=DZS (a, B)  • 
We  define  a  complex  1-D  array  TWAB(z)  =  w.  (We  use  the  prefix  T  in  a  FORTRAN 
name  to  refer  to  a  transformed  quantity,  as  the  caret  ~  is  used  with  an  algebraic 
name.)  For  a  pair  of  parameters  (a,B),  we  solve  the  complex  o.d.e.  (2)  for 
TWAB(z),  requiring  only  2-D  storage  in  (B,z)  on  the  right-hand  side.  The 
resulting  z-derivative  value  at  the  surface,  9w/9z,  is  first  evaluated  by 
second-order  one-s^ded  differences  on  TWAB(z),  and  then  stored  as  one  element 
in  DZS(a, B)  for  that  pair  (a, 3)  used  in  the  solution.  The  process  is  repeated 
for  other  (a,B)  pairs,  sweeping  out  the  soultion  w=TWAB (z) ,  and  again  requiring 
only  2-D  storage,  this  time  in  (a,B),  for  DZS  (a,B). 

The  solution  of  the  o.d.e.  (2)  can  be  achieved  directly  by  a  tridiagonal 

2 

solver  (e.g. ,  Appendix  A  of  Ref.  (13)),  presuming  that  0(Az  )  centered  differences 
are  used.  Higher-order  methods  such  as  O(Az^)  could  also  be  readily  incorporated, 
using  a  penta-diagonal  solver  for  O(Az^)  equations,  with  the  usual  ambiguity 
about  higher-order  boundary  conditions.  The  boundary  conditions  are  w=w=0 
at  z=-H  and  at  z=0. 

The  final  solution  is  obtained  by  reverse  transforming  the  doubly- 
transformed  variable  9w/9z=DZS (a, B) •  First,  the  3  transform  is  achieved 
directly  by  using  an  FFT  algorithm.  Second,  the  a  transform,  which  is  a 
Laplace  transform,  is  achieved  by  using  FFT  on  modified  data,  as  described 
by  Cooley  in  Ref.  (9).  This  reference  also  describes  a  method  for  determining 
the  arbitrary  parameter  of  the  Laplace  transform  from  considerations  of  accuracy 
and  the  maximum  x-value  desired.  The  solution  gives  9w/9z(x,y)  at  the  surface  z=0. 
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The  major  advantage  of  this  method  is  the  avoidance  of  the  eigenfunction 
solution  and  its  associated  stoi  age  problem.  The  only  disadvantage  is  the 
requirement  for  a  numerical  inversion  of  a  Laplace  transform.  This  can  be  a 
delicate  problem,  but  it  is  believed  that  the  methods  of  Ref.  (9)  will  permit 
accurate  inversion.  Other  accuracy  considerations  are  the  same  as  the  original 
method  proposed  in  Ref.  (1);  particularly,  we  note  the  error  associated  with 
the  assumption  of  periodicity  in  y,  required  by  the  use  of  a  discrete  Fourier 
transform  in  y.  It  is  expected  that  this  "aliasing"  error  will  be  small 
(Ref.  (12),  pg.  6). 

In  Ref.  (1)  (proposal,  pg.  25)  the  "exact"  boundary  condition  at  the 
free  surface  z=0  is  not  the  value  w(x,y  o)=0  used  above,  but 


w(x,y,o)  + 


(x,y,o)  =  0 
8x 


(4) 


It  is  stated  that  this  condition  "could  be  easily  incorporated,  if  needed" 
in  that  solution  method,  but  that  its  effect  is  expected  to  be  quite  small. 

In  order  to  incorporate  Eq.  (4)  in  the  presently  suggested  method,  an  over-all 
iteration  procedure  would  be  required.  The  entire  procedure  described  above, 
with  w(x,y,o)=0,  would  be  used  to  obtain  the  k=l  iterate  of  the  solution. 

Then  the  procedure  would  be  repeated,  using  as  a  boundary  condition,  for  the 
(k+1)  iterate, 


k+1 

w 


(x,y,o) 


U 82tj)k 
g  3x2 


(x,y,o) 


(5) 


Thus  the  incorporation  of  the  "exact"  boundary  condition  (4)  at  the  free 
surface  is  seen  to  be  possible  but  awkward  in  the  presently  suggested  method. 


Finally,  surface  velocities  are  solved  from  the  following  Poisson 
equation  (Ref.  (1),  proposal). 


.2,2  9z 

dx  dy 


(6) 
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The  velocity  potential  <j)  is  related  to  the  surface-wave  velocity  components 
by  u=3<j>/9x  and  v=9<j>/9y.  The  boundary  condition  at  x^  on  the  surface  is  set  by 
the  finite-difference  solution  to  the  near  wake.  (For  small  x^  or  large  11, 
the  turbulent  near  wake  will  not  have  "broken  through"  to  the  surface,  and 
<})=0  will  apply  at  x*^.)  At  y=0,  symmetry  prevails,  so  v=0.  For  the  other 
boundaries  in  y  and  x,  the  domain  must  be  large  enough  to  allow  the  use  of 
approximate  (computational)  boundary  conditions  on  (jj  or  9<y/9n,  such  as  the 
undisturbed  far-field  condition  (j)=0  or  the  less  restrictive  approximation 
9cf> / 3n— 0  (e.g.,  see  Ref.  (13),  Section  III-C).  The  solution  to  (3)  would  be 
achieved  with  any  of  several  fast  Poisson  solvers  (e.g.,  see  Ref.  (13), 

Section  III-B-1,8,9) .  Since  the  number  of  node  points  in  x  and  y  will  probably 
be  powers  of  2  (see  below)  and  the  boundary  is  rectangular,  a  good  choice  of 
methods  is  available;  since  FFT's  and  a  tridiagonal  solver  are  already  used 
in  the  solution  for  w,  Hockney's  method  (Ref.  (14))  is  a  natural  choice  for 
the  Poisson  solution.  Finally,  the  desired  surface-wave  velocity  components 
are  determined  by  simple  differentiation  of  the  solution  for  <}>,  with  u=9<J>/9x 
and  v=9<|>/9y. 

As  an  alternative  formulation,  it  was  suggested  in  Ref.  (1),  pg.  22, 
to  use  the  approximate  equation 


9v  _  9w 
9y  9x 


(7) 


at  the  surface.  This  simple  initial-value  problem  for  v(x,y)  can  be  solved  by 
marching  outward  in  y  from  the  symmetry  plane  y=0,  where  v(x,o)=0.  This  approxi¬ 
mation  is  very  easy  to  solve,  and  requires  no  further  approximations  of  computational 
boundary  conditions. 


4,2.3  Storage 

■k  *  * 

Since  the  complex  arrays  W  ,  dV  /dz,  R  (f3,z)  are  multiplied  by  complex 
coefficients  in  (2),  it  is  necessary  to  have  both  the  real  and  imaginary  parts 
in  core  at  the  same  time.  However,  the  full  complex  array  DZS(a,8)  may  be 
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stored  in  slow  core,  with  the  solution  values  buffered  out,  perhaps  a  line 
at  a  time,  through  a  singly-dimensioned  array  DABZS(B).  The  fast  core  storage 
required  at  one  time  would  be  the  3  complex  *-arrays,  and  the  two-singly- 
dimensioned  arrays  TWAB(z)  and  DABZS(B). 

The  storage  required  would,  of  course,  depend  on  the  resolution 
required.  For  256  points  in  the  £  (or  y)  direction  and  30  points  in  Z,  the 
^-arrays  require  6  *  256  x  30  =  46,080  words,  and  the  two  1-L)  arrays  require 
30  +  256  points  =  286,  for  a  total  of  46,366  words.  This  leaves  23,634  words 
of  fast  core  for  the  rest  of  the  program.  The  inversion  of  DZS  (ot,  3)  and  the 

Poisson  solution  for  <j)(x,y,o)  must  also  fit  into  fast  core,  but  these  can 

sequentially  share  the  locations  of  the  *-arrays.  For  64  points  in  the  a  (or  x) 

direction,  the  complex  array  DZS  (a, 8)  requires  2  x  256  x  64  =  32,768  words, 

which  is  less  than  the  ^-arrays  and  is  easily  accommodated.  (BLAND  COMMON 
would  of  course  be  used  to  utilize  the  array  storage  memory  for  compiling.) 

Going  entirely  to  the  500,000  words  of  slow  core  +  70,000  words  of 
fast  core,  we  could  have  tht  three  complex  ^-arrays  in  (8,z)  dimensioned 
(1024  x  65),  and  the  complex  DZS  (a, 3)  dimensioned  (64  x  1024),  giving  a  total 

of  6  x  1024  x  65  +  2  x  64  x  1024  =  399,360  +  131,072  =  530,432  words.  This 

would  leave  39,568  words  of  the  combined  fast  core/slow-core  memory  for  the 
rest  of  the  program. 

Other  combinations  will  occur,  but  it  is  desirable  that  the  number  of 
discrete  values  in  x  and  in  y  be  powers  of  2,  in  order  to  make  best  use  of  the 

speed  of  FFT  routines.  However,  a  single  factor  of  3  can  also  be  allowed  in 

some  versions  of  FFT,  without  completely  ruining  the  speed. 

4.2,4  Conclusion 

The  suggested  method  appears  to  have  advantages  over  that  described  in 
Ref.  (1).  No  eigenfunction  solutions  are  required,  and  only  two-dimensional 
arrays  need  he  stored.  A  numerical  inversion  of  the  Laplace  transform  is  required. 
The  suggested  method  is  an  alternate  to  the  approach  suggested  earlier,  and 
appears  to  be  preferable  to  it. 
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5.  FUTURE  WORK 

In  the  next  quarter  we  intend  to  continue  validation  of  the  code,  begin 
the  preparation  of  a  user's  manual,  and  to  concentrate  on  the  development  of  the 
far  field  part  of  the  code.  We  consider  the  completion  of  the  far  field 
capability  to  be  the  most  important  task  facing  us  at  the  present  time. 
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